﻿#include "LonosphereCorrection.h"
#include <cmath>
#include <iostream>
#include "TimeSys.h"
#include "CoordSys.h"
#include "literal.h" 


LonosphereCorrection::LonosphereCorrection() {
    // 初始化Klobuchar模型参数
    alpha[0] = 0.1397e-7;
    alpha[1] = 0.7451e-8;
    alpha[2] = -0.5960e-7;
    alpha[3] = -0.1192e-6;
    alpha[4] = 0.1192e-6;
    beta[0] = 0.1270e6;
    beta[1] = 0.1966e6;
    beta[2] = -0.6554e5;
    beta[3] = -0.2621e6;
    beta[4] = 0.1966e6;
}


LonosphereCorrection::~LonosphereCorrection() {
}

double LonosphereCorrection::calculateIonosphereDelayKlobuchar(const BaseT& time, const XYZ& receiverPos, const XYZ& satellitePos) {
    // 将接收机的XYZ坐标转换为BLH坐标
    BLH receiverBLH = CoordSys().XYZ2BLH(receiverPos);
    // 调用另一个重载函数进行计算
    return calculateIonosphereDelayKlobuchar(time, receiverBLH, satellitePos);
}

double LonosphereCorrection::calculateIonosphereDelayKlobuchar(const BaseT& time, const BLH& receiverPos, const XYZ& satellitePos) {
    // 使用 literal.h 中定义的 OMEGA_PZ90
    double omega_e = OMEGA_PZ90;

    double delay = 0.0;

    // 将接收机的BLH坐标转换为XYZ坐标
    XYZ receiverXYZ = CoordSys().BLH2XYZ(receiverPos);
    // 计算卫星和接收机之间的差值向量
    XYZ diff = { satellitePos.x - receiverXYZ.x, satellitePos.y - receiverXYZ.y, satellitePos.z - receiverXYZ.z };
    // 计算卫星和接收机之间的距离
    double distance = sqrt(diff.x * diff.x + diff.y * diff.y + diff.z * diff.z);
    // 计算卫星的仰角
    double elevation = asin(diff.z / distance);

    // 计算时间参数
    double t = static_cast<double>(time.sec) + time.sesec;
    double t0 = fmod(t, 7200);
    if (t0 < 0) t0 += 7200;

    // 计算纬度和经度的弧度值
    double psi = receiverPos.L * D2R;
    double lambda = receiverPos.B * D2R;

    // 计算倾斜因子
    double F = 1.0 - 0.00266 * cos(2.0 * lambda) - 0.00028 * receiverPos.H / 1000.0;

    // 计算Klobuchar模型中的A和B参数
    double A = alpha[0] + alpha[1] * cos(2.0 * psi) + alpha[2] * sin(2.0 * psi) + alpha[3] * cos(4.0 * psi) + alpha[4] * sin(4.0 * psi);
    double B = beta[0] + beta[1] * cos(2.0 * psi) + beta[2] * sin(2.0 * psi) + beta[3] * cos(4.0 * psi) + beta[4] * sin(4.0 * psi);

    // 计算总电离层延迟
    double T = A + B * cos(omega_e * t0);

    // 限制T的范围在0-100纳秒之间
    if (T < 0) T = 0;
    if (T > 100e-9) T = 100e-9;
    
    delay = T * F;
    
    // 将结果从秒转换为米(光速约为2.99792458e8 m/s)
    return delay * 2.99792458e8;
}

void LonosphereCorrection::test() {
    // 定义测试时间
    BaseT time = { 1491534395, 0.0 };
    // 定义测试接收机位置
    BLH receiverPos = { 34.808610734, 113.357448942, 120.2032 };
    // 定义测试卫星位置
    XYZ satellitePos = { -32348377, 27042006, 509548 };

    // 计算电离层延迟误差
    double delay = calculateIonosphereDelayKlobuchar(time, receiverPos, satellitePos);
    // 输出电离层延迟误差改正值
    std::cout << "电离层延迟误差改正值: " << delay << " 米" << std::endl;
}
